The field of tomography involves obtaining a set of signal measurements (such as X-ray attenuation or traveltimes of acoustic or electromagnetic waves) that traverse a physical domain and then employing an inverse method to estimate spatially distributed material parameters (such as absorption or propagation velocity). The results are maps of spatially distributed estimates of the material properies, which are typically displayed in graphical format. Tomographic methods were first developed for biomedical imaging applications, such as disclosed in U.S. Pat. No. 3,778,614 to Hounsfield. These methods were subsequently adapted for geophysical applications using cross-borehole measurements as described in the article by Dines et al., entitled "Computerized Geophysical Tomography," Proc. of the IEEE, 67(7):1065-1073, 1979.
Although there have been many improvements in tomographic methods since they were first introduced, shortcomings remain. In particular, there remains a need for accurate and computationally efficient algorithms to invert possibly sparse three-dimensional signals, as well as a means for assessing the spatially-varying accuracy of the inversion. While many so-called three-dimensional tomographic methods exist, the vast majority of these methods actually invert the data in two-dimensional planes, and then construct the three-dimensional image from these two-dimensional slices. In the publication by Miller et al., entitled "3-D maximum a posteriori estimation for single photon emission computed tomography on massively-parallel computers," IEEE trans. on Med. Imaging, 12(3):560-565, 1993, it is shown that full three-dimensional single photon emission inversion is significantly more accurate than stacked two-dimensional inversions. However, the resulting method is so computationally intensive that it requires a massively parallel implementation. The method described therein employs an iteratively applied back projection method, which requires dense measurements. Similarly, in the publication by Alumbaugh et al., entitled "Three-dimensional massively parallel electromagnetic inversion--II. Analysis of a crosswell electromagnetic experiment," Geophys. J. Int., 128:355-363, 1997, a massively parallel full three-dimensional geophysical electromagnetic inversion using the iterative conjugate gradient method is disclosed.
There are trade-offs inherent in specifying the degree of spatial discretization in tomography. Very finely discretized models require more computational resources and are typically implemented using relatively simple iterative update mechanisms which may require additional arbitrary interpolation, smoothing, and regularization. Furthermore, unless data density and quality are similarly high, fine discretization can lead to instabilities in the inversion. On the other hand, coarsely discretized models, which are more computationally efficient and stable, may not provide adequate resolution of the heterogeneity present and detectable in the signal, and errors in the geometry of the parameterization may lead to errors in the estimation.
Unlike typical tomographic methods, which estimate model parameter values only, stochastic inversion methods model discretized parameters with stochastic random variables, in which each discrete parameter is modeled by its mean value, its variance, and its covariance with all other parameters in the domain model. Perhaps the most statistically sophisticated stochastic inversion algorithm is the extended Kalman filter, which has been widely and successfully applied to problems in control systems engineering. See, e.g., Gelb, A. (Ed.), Applied Optimal Estimation, The M.I.T. Press, Cambridge, Mass., 1974. This Bayesian filter uses prior estimates of the model values and error covariance of the parameters and the states, and the state-parameter cross covariance. The filter updates the model values and covariance in accordance with measurements so as to minimize the estimation error variance of the estimated values, by explicitly estimating the spatially-varying measurement error covariance and system noise covariance of the forward simulator. The updates are optimal for the linearized system. A pseudo-code description of the extended Kalman filter is shown below.
Time Update:
x=.PHI.x+f (i.e., predict states x from parameters y using a computer programmed forward simulator) ##EQU1## EQU P.sup.xy =JP.sup.y +.PHI.P.sup.xy EQU Px=JP.sup.y J.sup.T +Q+.PHI.P.sup.x .PHI..sup.T +J(P.sup.xy).sup.T .PHI..sup.T +.PHI.P.sup.xy J.sup.T
Measurement Update: EQU K.sup.x =P.sup.x H.sup.T (R+HP.sup.x H.sup.T).sup.-1 EQU K.sup.y =(P.sup.xy).sup.T H.sup.T (R+HP.sup.x H.sup.T).sup.-1 EQU x=x+K.sup.x (z-Hx) EQU y=y+K.sup.y (z-Hx) EQU P.sup.x =P.sup.x -K.sup.x HP.sup.x EQU P.sup.y =P.sup.y -K.sup.y HP.sup.xy EQU P.sup.xy =P.sup.xy -K.sup.x HP.sup.xy
wherein
= denotes variable assignment in a computer program; PA0 .PHI. is the transition matrix (a function of y) for predicting x; PA0 x is a vector of state variables PA0 y is a vector of parameter variables; PA0 J is the Jacobian sensitivity matrix; PA0 f is a vector of forcing conditions; PA0 P.sup.x is the state covariance matrix; PA0 P.sup.y is the parameter covariance matrix; PA0 P.sup.xy is the state-parameter cross-covariance matrix; PA0 Q is the covariance of the system noise of the simulator; PA0 R is the covariance of the measurement error; PA0 H is the observation matrix; PA0 z is a vector of measurements; PA0 K.sup.x is the state gain matrix; PA0 K.sup.y is the parameter gain matrix; PA0 superscript T denotes matrix transpose; and PA0 superscript -1 denotes matrix inversion.
The covariance matrices provide implicit smoothing, interpolation, and regularization. One or more measurement updates may be performed during a given timestep. The recursive updates of the parameter and state values and covariances enable sequential incorporation of time-series data measurements, as well as sequential assimilation of subsets of batches of stationary measurements (to limit the size of the matrix inversions required) and/or iteration upon the same data to help reduce errors due to non-linearities. The resulting covariance estimates provide a means of assessing the accuracy of the parameter estimates. The computational complexity of the extended Kalman filter has limited its application to the estimation of relatively few parameters, and it is believed this method has not previously been used in tomography.
Iterative co-kriging is described in the publication by Yeh et al. entitled "An iterative stochastic inverse method: Conditional effective transmissivity and hydraulic head fields", Water Resources Research, 32(1):85-92, 1996. Iterative co-kriging is similar to, but less general than, the extended Kalman filter. In iterative co-kriging, measurement error and system noise are not explicitly modeled. Instead, a regularization term is added, and parameter (but not state) values and covariance can be conditioned on the measurements. Sequential co-kriging, as described in the publication by Harvey et al., entitled "Mapping hydraulic conductivity: Sequential conditioning with measurements of solute arrival time, hydraulic head, and local conductivity," Water Resources Research, 31(7):1615-1626, 1995, is similar to iterative co-kriging, but designed to be used to sequentially incorporate distinct forms of data. Co-kriging methods have been applied extensively to the two-dimensional groundwater inverse problem, but not to very large problems as are typical in three-dimensional tomography.
Cluster analysis is the formal study of algorithms and methods for grouping, or classifying, objects (see Jain et al., Algorithms for Clustering Data, Prentice-Hall, Englewood Cliffs, N.J., 1988). In the well-known k-means clustering algorithm, data are organized into one of k clusters, each data element being placed in the cluster whose cluster centroid (i.e., the mean of all data elements in the cluster) is closest to the value of the element. The number of clusters can be predetermined or can be based on the characteristics of the data itself. For example, divisive, partitional clustering algorithms cluster the data into as few clusters as possible, such that the data elements within each cluster are more similar to each other than to data elements in other clusters. Clusters are added in as necessary to satisfy a particular clustering criterion. Several imaging methods, such as that disclosed in U.S. Pat. No. 4,751,643 to Lorensen et al., have been devised using cluster analysis methods to post-process three-dimensional reconstructed tomographic images to identify distinct (e.g., anatomical) structures of a larger scale than the spatial discretization of the imaged parameter. The article by Hyndman et al., entitled "Traveltime inversion for the geometry of aquifer lithologies," Geophysics, 61(6):1728-1737, 1996, discloses a method wherein two-dimensional cross-borehole traveltime tomography is alternated with cluster analysis in an iterative fashion to determine aquifer lithology. The method involves constraining the assignment of pixel-based parameter values to one of a small number of populations based on the histogram of slowness residuals, thereby increasing stability and convergence properties (but also the computational complexity) of the algorithm. Two-dimensional images are subsequently interpolated into three dimensions. In these non-stochastic methods, parameter covariances are not modeled. Therefore, when parameters within clusters are merged, there is no need to merge their covariances. These structure determination methods do not employ clustering as a means for estimating the parameterization for subsequent estimation. Therefore, the cluster analysis increases, rather than decreases, the computational burden of these imaging technologies.
Statistically homogeneous two-dimensional random fields can be merged together using a technique known as random field averaging. See, e.g., Vanmarcke, E. Random Fields. Analysis and Synthesis, The M.I.T. Press, Cambridge, Mass., 1983. The variance of the averaged field is always less than or equal to the area-weighted average of the variance of the component fields, even when the estimates of the means of the component fields are different. If the true means of the component fields are different (i.e., they are statistically heterogeneous), then random field averaging will underestimate the true variance of the merged field. This method is therefore of limited usefulness for "upscaling" heterogeneous domains modeled by heterogeneous stochastic random variables.